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Abstract 

In this paper the multi-dimensional random walk models governed 
by distributed fractional order differential equations and multi-term frac- 
tional order differential equations are constructed. The scaling limits of 
these random walks to a diffusion process in the sense of distributions is 
proved. Simulations based upon multi-term fractional order differential 
equations are performed. 
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1 Introduction 

1.1 Motivation. 

In this paper we study simulation models based on distributed order differential 
equations, which we will call DODE simulations. This type of simulation reflects 
the rich structure of diffusion media, in which a several diffusion modes are pos- 
sible. Diffusion processes with complex and changing modes are ubiquitous in 
nature (see, [51 151 1221 1261 13U| and references therein). One of the motivations for 
conducting DODE simulations is to model the movement of proteins on the cell 
membrane. Numerous experiments |1(JI II 11 [TBI 1241 125| show that macromolecule 
movement through the cell membrane is distinct from Brownian motion. Saxton 
and Jacobson [25] noted that practically all experimental results show apparent 
transitions among modes of motion. 
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The governing equation, which we take as a basis for our simulation models, 
in general form, is distributed space fractional order differential equation 

D%u(t,x)= [ a{a)D%u(t,x)da, t> 0, x € M N , (1) 
Jo 

where < j3 < 1 , D* is the Caputo fractional order derivative ^] , Dq = 
(—A) a is the space fractional order (pseudo-differential) operator with the sym- 
bol |£| a . Note that Dq can be written in the form of hypersingular integral as 
well [23 • The function a(a) is a positive integrable function (or positively de- 
fined distribution). Depending on a(a), Q may become a multi-term fractional 
order differential equation, which can possibly describe the existence of a finite 
number of diffusion regimes. Although, the distributed order differential oper- 
ators were first mentioned by @] EJ in the f 960s, the intensive study of models 
based on the distributed order differential equations has been started recently 
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The present report is organized as follows. In Section 2, we briefly recall 
the theoretic platform of the construction of the DODE simulation models an- 
nounced in • In Section 3 we analyze the difference schemes associated with 
the DODE models, and in Sections 4 and 5 we construct random walk models 
and simulations based on the transition probabilities introduced in the previous 
sections. 

1.2 Notation. 

In this paper, 1R N is the iV-dimensional Euclidean space with coordinates x = 
(xi, ...,xn) while Z N is the ^-dimensional integer-valued lattice with the lattice 
nodes being given by the multi-index notation j = (ji, jw). The letters i, j 
and k will be exclusively used for the multi-indexing of lattice nodes. We denote 
by Xj — (hj lt hj N ),j £ Z N , the nodes of the uniform ^.-lattice which is 
defined as (hZ) N with h being the distance between any two lattice nodes. We 
introduce a spatial grid {xj = jh,j £ Z N }, with h > and a temporal grid 
{t n = nr, n = 0, 1, 2, ...} with a fixed stepsize r > 0. Furthermore, let u™ denote 
the discretization of the function u(t, x) on the spatial and temporal grid at 
x = Xj and t = t n , i.e «" = u(t n , Xj). 

2 Markovian random walks associated with the 
DODE 

2.1 Particle jumps. 

Assume X to be a N-dimensional random vector [3T| whose values range in 
Z N . Let a sequence of random vectors Xi, X2, ... also be N-dimensional inde- 
pendent identically distributed random vectors, all having the same probability 
distribution. Consider the sequence of random vectors 

S n = hX-i + hX-2 + ... + /iX„, n — 1, 2, ... 
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taking S = = (0, . . . , 0) G for convenience. We interpret Xi, X 2 , as a 
sequence of particle jumps starting time t 
takes a jump h~K n from S„_i to S„. 
particle being at location Xj at time t. 
S n+ i = S n + /iX n+ i, we have 



to = 0. At time t — t n , the particle 
u(t n ,Xj) is the probability of a 
and, taking into account the recursion 



If u- 



k£Z N 



3 EZ N , n = 0, 1, 



fc) 



(2) 



where the coefficients pk, k € Z N are called the transition probabilities. The 
convergence of the sequence S„ when n — > oo means convergence of the dis- 
crete probability law (probability mass function) (u") jeZ jv, properly rescaled 
as explained below, to the probability law with a density u(t, x) in the sense of 
distributions (in law). This is equivalent to the locally uniform convergence of 
the corresponding characteristic functions (see for details PU). This idea is used 
in |271 123) to prove the convergence of the sequence of characteristic functions 
of the corresponding random walks to the fundamental solution of distributed 
order diffusion equations. 

2.2 Markovian transition probabilities. 

Let the transition probabilities in Eq.@ take the form 



where 



Qk{a, h) 





Pk 


f 


a{a)b{a) 


10 


_\k\ N + a h a _ 



da, and b(a) 



[r(i + §)] 3 



7T 2 2 



The transition probability pq can then be defined as 

Po = 1 - E Pfc = 1 - T Qo( a > h )> 

k^O 



where 



(a, h) = ^2 1k(a, h) = ^2 j 

k^O k^Q J ° 



f 


a(a)b(a) 


10 


_\k\ N + a h a _ 



da, 



(3) 



(4) 



(5) 



(6) 



Assuming that the condition < Tqo(a, h) < 1 is fulfilled, the transition prob- 
abilities then satisfy the following properties: 

1- Pfc = l; 

kez N 



2. Pk > 0,k e Z N . 
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Note that the non-negativity condition 2 in property 2 is linked with the Riemann 
zeta-function. Indeed, introduce the function 



1 00 M 

R(a ),y ' ^J^, o<«<2, 



(7) 



where M m — J^|fe|=m ^ n one-dimensional case IZ(a) = 2£(1 + a), where 
£(z) is the Riemann zeta-function. Then the nonnegativity condition < po < 1 
can be rewritten as 



rq (a, h)=T I 

Jo . 



a(a)b(a)lZ(a) 



da < 1. 



(8) 



It follows from this condition that h — > yields t — * 0. This, in turn, yields 
t/r —> oo for any finite t. 

Theorem 1 Let X be a random vector with the transition probabilities pk — 
P(X = Xk),k € Z N , defined in Eq. J3J) and i?g. (0) which satisfy properties 1 
and 2. Then the sequence of random vectors S„ = ftXj + ... + /iX n , converges 
as n — > co m tow to t/ie random vector whose probability density function is 
the fundamental solution of the distributed space fractional order differential 
equation 0) with /3 = 1. 

Note, for the simulations used in this paper, it is important to use the multi- 
term analog of this theorem. Assuming that 



M 

a(a) = 2J a m S(a - a m ), < cti < ■ ■ ■ < a M < 2, 

m— 1 

with positive constants a m , we get a multiterm DODE 

M 

D^u(t,x) = a m D% m u{t,x), t > 0, x e M N . 

rn — 1 



(9) 



(10) 



Also note that the coefficients flfe(a, h) in Eq.QJ and Eq.@ become multi-term 
as well: 



M 



qk(a,h) = Y 



a m b(a m ) 
\k\ N+a ™h a « 



Theorem 2 Lei the transition probabilities pk = P(X = Xfc),fc G Z N , 0/ i/ie 
random vector X &e given as follows: 



Pk = rq k (a, h) and p = 1 - Tq (a, h) 



(11) 



2 This condition is equivalent to the stability condition of finite-difference schemes giving 
the usual stability condition if a(a) = &{a — 2). 
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where a(a) — 5Zm=i a m5{a — a m ). Assume 

Ea m b{a m )1Z{a m ) 
; < 1- 

m— 1 

Then the sequence of random vectors S n = ZiXi + ... + /iX„, converges as n — > oo 
in law to the random vector whose probability density function is the fundamental 
solution of the multiterm fractional order differential equation hl(J\) with (3 = 1 . 

Remark: 

As we noted above these results were announced in ]2ty . The more general 
case of these theorems corresponding to a fractional (3 £ (0, 1) can be 
obtained introducing a positive waiting time distribution and corresponding 
iid random variables [13. 20L We do not describe this case in this paper. 
We note only that the general case is studied by applying a general finite- 
difference approach and that this general difference scheme is stable under 
some condition and has a unique solution. 



3 Generalized Transition Probabilities for the 
DODE 

The set of grid points in used to update u at time t = t n+ \ = (n + 1)t 
is called the stencil. In this section, we start from stating the values of the 
transition probabilities associated with the stencil for the discretization of the 
particular space-time-fractional differential equation, 

D%u(t, x) = D%u(t, x), t > 0, x 6 M N , < (3 < 1, < a < 2, (12) 

and then generalize it to distributed order differential equations. 



3.1 Discretization of the time- fractional derivative. 



Using the Caputo time- fractional derivative 0, the left-hand-side of (|12|l be- 
comes 



du(s, x) 



ds 



ds 



(t- s y 



<[3<1. 



(13) 



Note that when = 1, D%u(t,x) = du/dt. When < [3 < 1, we will use the 
following discretization (see rT7] for the derivation): 



1 



r(i - 13) 



E 



uAt n+1 - s) 



ds 



r 1 - E 



(14) 



711 — 1 
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where 



In 



1 



,1-/3 



TO 1 ^, TO 



0,1,. 



7m-l - 7m, "2 



.,71 



and i/ = r(2 — j3). The formulas for the coefficients c m and 7 m and the scalar 
v that were used in i|14|) . which were based upon the Caputo time- fractional 
derivative, easily generalize to other definitions of the time-fractional derivative. 
For example, in the case of the Grunwald-Lctnikov time-fractional derivative, 
v = 1 and 7 m and c m are re-defined as the following 



rn 



1. 



7» 



1 — y^Cj, m = 0, 

i=l 



For simplicity of notation, we will now set 



w 



In 

Cn+i — 1 , ^ 1 , • 



and, as a result, (|14|) can be rewritten as 

1 




w m u, 



Note that for /? = 1, z/ = T(2 - /3) = 1 and w 



(15) 



w n _i = with w r , 



In this case, (|14l) reduces to the standard forward-time discretization for du/dt: 



du 



,n+l 



3.2 Discretization of the space- fractional derivative. 

Just as the discretization for the time-fractional derivative assumes a simple 
form when = 1, the discretization for the space- fractional derivative, based 
upon centered differences, assumes a simple form when a = 2. For example, 
when a — 2 and the N = 2, 

D >] = ^l~h +«WiA+D +U U^-D 

In TV-dimensions, the stencil consists of j = (ji, ■ ■ ■ ,Jn) an d hs nearest 2iV 
neighbors with each nearest neighbor being h units away from j. When a = 
{ai, . . . , cum} 7^ 2, the space-fractional derivative is given by |29|: 

Dgv% « -ffl,(a, fc)u7 + ^ % (a, /i)uj_ fe (16) 

fe#0 

where the coefficients qo(a, h) and <?fc(o!, /i) are defined in Q and © using the 
multiterm definition for a(a). The geometric consequence of changing a from 
a = 2 to a = {ai, . . . , au} ^ 2 is that the stencil gets enlarged from 2N + 1 
grid points to all of the lattice points in . 
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3.3 Construction of the explicit finite difference scheme. 

Setting the discretizations for the time and space-fractional derivatives equal to 
each other in l|14|l and (|16[) . we get 

^ w m uf \ = -q Q (a,h)u% + ^q k (a,h)u]_ k . (17) 
Solving for itf , the following explicit finite-difference scheme is constructed: 




n-l 

m=0 keZ N 



a': ' = > tr,„<iT+ > Pkuf_ k , (18) 



where 

p k = vT^Qkia, h), k ^0 and p = w n - ur l3 q (a 1 h). 

When [3 — 1, the coefficients pk arc equivalent to the transition probabilities pk 
in |XJ) . Furthermore, since all the transition probabilities are non- negative and 
taking into account that w n = c\ = 2 — 2 1_/3 and ^ = Gamma(2 — (3), we have 
an upper bound for the stepsize r: 

/ 2-2 1 -' 3 x 1/13 

p Q >0 < r < 

VT(2 - f3)q {a,h) 

The update in (|18(l is determined by Markovian contributions (those 

values of u at time t = t n ) and non-Markovian contributions (those values of u at 
times £ = {to, t\, . . . , i n ~i}). The order of the time fractional derivative (3 deter- 
mines the effect that the non-Markovian transition probabilities (wq, . . . , w n -\) 
has on u" +1 . This effect can be measured by examining the sum of all of the 
transition probabilities in (|18(l : 



n-l 

w m + ^2 Pk = !> 

m=0 k£Z N 



n-l 

= 1 - u>„ 



kez N 



(19) 

Pk = w n . 



Recall that when /3 = 1, w n = 1 and wq = ■ ■ ■ = m,,-i =0. In this case, the 
first term in (TH?)l vanishes and p = 1 — Tq (a, h). 

When < [3 < 1, the values of u™ associated with t S {to, ■ • ■ i^n-i} are 
weighted by the coefficients {wq, w%, . . . , u> n _x}. Figure plots w m for m = 
0,1,..., n where n — 100 and (3 — 0.9. It is well-known that the sequence 
{ w m}m=i are monotone increasing [8], i.e. Wi < u>2 < ■■■ < w n ^i < w n . 
However, it is not true wq < w%. In fact, in Figure^ wg < wq < w%. Hence, 
the contribution of w° to wj 01 is quite large relative to the other intermediate 
values of u™. We will see later on that this will have important consequences in 
non-Markovian random walk numerical simulations. 
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Figure 1: The weight w m associated with the density it™ is plotted as a function of 
m for both the Caputo and Griinwald-Letnikov (G.L.) time-fractional derivatives and 
j3 — 0.9. The lower dotted horizontal line corresponds to the value of wo ~ 0.005 while 
the upper two dotted lines correspond to w n = ci for both the Griinwald-Letnikov 
(iwioo = 0.8) and Caputo derivatives (wioo ~ 0.851). 

4 Monte Carlo Protocol for the Random Walk 

4.1 General Framework. 

The random walk model corresponding to the governing equation in (|12fl uses 
the non-Markovian transition probabilities, w m the the Markovian transition 
probabilities pk to assign where in the Zj^ lattice a particle will jump to. This 
jump can be based upon a partitioning of the unit interval V — [0, 1) into two 
disjoint subintervals V\ and Vi such that V = V\ U Vi where V\ = [0, 1 — w n ) 
and V% = [1 - w n , 1). 

We will use a two-dimensional walk for illustration purposes. The random 
walk process begins by generating a uniformly distributed random number r 
in the unit interval and observing what subinterval (Vi or Vi) it falls into. If 
r e Pi = [0, 1 — w n ), then the particle will do a non-Markovian jump, i.e. 
the jump will be determined by transition probabilities w m , m = 0, . . . , n — 1. 
Otherwise, if r 6 V = [1 — w n , 1), then the particle will undergo a Markovian 
jump, i.e. the jump will be determined by transition probabilities p^. In effect, 
the random walk interpretation presented here is a two-dimensional extension 
of the one-dimensional random walk interpretation given in |14|. 

4.2 Non-Markovian Jumps. 

If < (3 < 1 and r S V\, then the jump that the particle takes will be determined 
by w m ,m = 0, . . . ,n— 1. Let A = {Ao , At , . . . , A n — i } be an n-element set such 
that Ai — Wi, i — 0, . . . , n — 1. Furthermore, let the interval V\ be refined in 
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the following way: 

V 1 = [B ,B 1 ,...,B n ), 

such that Bq — and Bj = J2iZo ^i, j = 1, . . . , n. If r S [Bo, Si) = [0, wo), 
then the position of the particle at t = t n+ i is given by S n +i = So (the origin). 
Otherwise, if r £ [Bj—i, Bj), j — 1, . . . , n, then the particle will jump back to 
the position that it had visited at time t = tj, i.e. S n +i = Sj. 

4.3 Markovian Jumps when a = 2. 

If r 6 7*2 = [1 — io n , 1) and a = 2 then the jump will only be to adjacent lattice 
grid points. Let Vn be partitioned in the following manner: 

V 1 = [B ,B 1 ,...,B 5 ) 

where B = 1 — w n and Bj = B + J2i=o C? = 1, •••>5). Here, *4 = 
{A), *4.2, -^3, Ai} where Ao = w n - Ar\ and A% = r\ = vt® jh a , i = 1,2, 3, 4. 
If r £ [Z?o,Bi), then the particle remains at the current position, otherwise if 
r G {[#1,62), [#2,63), [^3,^4), [£>4,£>5)} then the particle will move left, right, 
up or down, respectively, one lattice position. 

4.4 Markovian Jumps when a = {ocx, . . . , %} 7^ 2. 

If r G 7^2 = [1 — 1) and a = {ai, . . . , aw} 7^ 2, then the jump will be 
determined by an infinite partition refinement of T>2- Let 

A = {A ,Ai,...}; Vi = [B Q ,B x ,...) 

such that Bo = 1 — w n and Bj = Bq + J2i=o "4* (j = 1, • • ■)■ I n this case, the 
set A consists of all of the transition probabilities pk, k G Z 2 , with Ao = po- 
If r G [Bo,Bi) = [1 — w n , (1 — w n ) + po), then the particle will remain at 
the current position. Otherwise, if r G [B s ,B s+ i), then there exists a unique 
k = (fci, £2) G associated with s G W such that the particle will jump from 
S„ to S n+1 = S„ + (kih, k 2 h). 

5 Simulations 

Our motivation of the numerical simulations presented here is to see how DODE 
simulations of biomolecular motion of particles on a cell surface differ from those 
based upon classical Brownian motion. Although the DODE random walk mod- 
els are described theoretically for multivariate case in iV-dimensions, neverthe- 
less all our simulations are conducted in the two dimensional case since we are 
interested in the diffusion of proteins on a cell membrane surface, which can be 
locally approximated by a two-dimensional membrane sheet. In |15|. simulated 
particle motion is based upon the classical Brownian motion scenario (where 
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a = 2 and (3 = 1) in which the particle is confined within cytoskeletal bar- 
riers (see Figure 01. In these single particle tracking studies, particle appears 
to be spatially and temporarily confined within transient confinement zones. 
Although the barriers are never directly observed, it is postulated that the cy- 
toskeletal barriers are the reason for the transient spatial confinement of particle. 
In principle, DODE simulations provide an alternative explanation for the ob- 
served trajectories in single particle tracking studies that does not necessarily 
require the existence of cytoskeletal barriers to explain transient confinement. 

In |TS], the authors use the mean-squared-displacement formula 4ar = h 2 
in which the parameters a (the diffusion coefficient) , r (the timestep) and h 
(the lattice width), respectively, are given using the following values: h = 6 
nanometers and r = 1/is (microseconds, or r = 1CP 6 seconds). Since the mean- 
squared displacement formula implicitly assumes that 

p = l-4a^ = l-4al = 0, 

the diffusion coefficient is then computed as a — h 2 /(4r) = 9 x 10~ 12 m 2 /s. 
To facilitate a comparison of our DODE simulations with the simulations of 
|15l lift) , we will also use the same diffusion coefficient (a\ = ■ ■ ■ = « M = a = 
9 x 10~ 12 m 2 /s) and the same lattice width (h = 6 nanometers). Using the fact 
that the transition probabilities sum to 1, 

n-1 

1 = ^ W m + ^Pk = {1 - w n ) +p + VT P q a {a, h) 

m=0 k 

we can now solve for r in terms of a, (3 and p , 



\uq (a,h)J \T(2 - (3)q (a, h) 



As in US]) we set Po = 0. However, due to the dependence of t on a and /3, the 
relative size of the timestep (from r = 10~ 6 s in the case of a = 2 and (3=1) 
will change as a and (3 vary. Instead if fixing the simulations to have the same 
stepsize r, we will fix the duration of the overall walk to be the same, Let T 
denote the overall duration of the random walk simulation. In all of our DODE 
simulations, T is set to T = ^ seconds. This is equivalent to 1 frame at video 
rate where video rate is measured as 30 frames per second. All simulations were 
performed in M ATL AB [T§] . 

Figure |3 shows various Markovian DODE simulations ((3=1) across various 
values of a. The left, middle and right plots in the top row show DODE simula- 
tions for a = {2}, a = {1.5} and a = {1.5, 2}, respectively. The first two DODE 
simulations are actually monofractal DODE simulations with M = 1 while the 
last one (a = {1.5, 2}) is a multi-fractal case with M = 2. The large white dots 
indicate the first and last positions of the random walk and the starting posi- 
tion is always the origin (0, 0). It is clear that for these DODE simulations with 
a =/= {2} that the particle travels much longer distances since the probability 
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Figure 2: This random walk simulation depicts classical Brownian motion confined 
to rectangular cytosketetal barriers. The parameters used in this simulation are as 
follows: h — 6 nanometers, r = 10 _6 s and a = 9 x 10~ 12 m 2 /s. The barriers are 
spaced out every 66 nanometers and the the probability of escape is p — 0.01 when a 
particle encounters a barrier. 

of jumping to faraway lattice sites is greater than what would be expected for 
a = 2. 

Figure 0] shows various non-Markovian DODE simulations ((3 = 0.999) us- 
ing the same values of a as in Figure EH The bottom plot in both Figures 
[3] and ^ show the plots on top row superimposed on one graph. The dark 
shaded lines correspond to Markovian jumps (r G V\) while the white lines 
correspond to non-Markovian jumps (r G V%). The frequency of the non- 
Markovian jumps are given by the size of the V\ interval. For (3 = 0.999, 
V\ « [0, 1 — w n ) — [0,0.00069339). Hence, the probability at every timestep 
of doing a non-Markovian jump is 0.00069339. The bottom plot in Figure 0] 
shows the superposition all three non-Markovian DODE simulations on the same 
graph. 

For Figure El we have non-Markovian DODE simulations for a fixed set of 
a values (a = {0.8, 1.3, 1.8}) with (3 varying. The left, middle and right plots 
correspond to (3 = 0.999, (3 = 0.99 and (3 = 0.9, respectively. The probability 
of taking a non-Markovian per timestep for these graphs is 0.00069339 (left), 
0.0070 (middle) and 0.0718 (right). For example, roughly 7% of all jumps for the 
right subplot on the top row are non-Markovian jumps. The effect of decreasing 
(3 is clear: the overall distances that the particle traverses is decreased since 
motion is constrained by jumps to previously visited positions. 

The average jump sizes associated with Figures[310]and[S]are shown in Table 
2] The numbers in the brackets before the colon correspond to the (a, f3) pair 
used in the DODE simulation while the number after the colon corresponds 
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to the average jump size. For the non-Markovian walks, the average jump 
length is larger when, for a fixed set of a values, (3 is decreased from 1. This 
is a consequence of the non-Markovian nature of the random walks for < 
(3 < 1. Since the particle is allowed to jump back to any previously visited 
position, the jump size can be quite large if the previously visited position 
was spatially remote from the particle's current position (see Figure |SJ . In 
particular, in Figure^ the probability of the particle to jump back to the origin 
is disproportionately larger than for other previously visited sites. In Figures 0] 
and |S1 one can observe evidence of this phenomenon. 



Table 1: This table reports the average jump size (after the colon) for all of 
the DODE simulations in Figures E] and The numbers before the colon 
indicates values of the (a, /3)-pair used in the DODE simulation. 





Left Plot 


Middle Plot 


Right Plot 


Figure 


(2,1): 6.0000 


(1.5,1): 10.9770 


({.5,2},1): 7.3320 


Figure 


(2,0.999): 6.0038 


(1.5,0.999): 11.0707 


({1.5, 2}, 0.999): 7.3593 


Figure |3 


({0.8,1. 3,1.8}, 0.999): 17.0328 


({0.8, 1.3, 1.8}, 0.99): 17.1663 


({0.8, 1.3, 1.8}, 0.9): 19.8946 



6 Conclusion 

Qualitatively, the DODE simulations provide a richer repertoire of motion, com- 
pared to monofractal walks when M = 1. Macroscopically, the DODE trajecto- 
ries tend to cluster together more often than the monofractal walks. The clus- 
tering is even more pronounced when the motion is non-Markovian due to the 
memory the particle has for previously visited positions. Moreover, one does not 
have to hypothesize the existence of barriers to explain why a particle appears 
trapped in a transient confinement zone or hops large distances. The clustering 
of trajectories and large jumps are a natural consequence of the DODE random 
walk model. However, when the motion is non-Markovian, the particle has a 
strong propensity to jump back to the origin, a consequence of the dispropor- 
tionately large weight wq associated with u®. While jumping back to previously 
visited "compartments" is observed for experimentally observed single particle 
tracking data one does not experimentally observe molecules jumping back 
from its current position to the starting point. Nonetheless, the DODE random 
walk models closely resemble the data from single particle tracking experiments 
of molecules moving on cell membranes |1 51 116) . This is not surprising since 
the motion of biomolecules on the cell surface occurs in a very heterogeneous 
environment. 
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Figure 3: The first three subplots in the top row correspond to Markovian DODE simula- 
tions (/3 = 1) with different values of a: a = 2, a = 1.5 and a = {1.5,2} for the left, middle 
and right plots. The bottom plot superimposes all of the top three simulations on one graph. 
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Figure 4: The first three subplots in the top row correspond to non-Markovian DODE 
simulations (/3 = 0.999) with different values of a: a = 2, a = 1.5 and a = {1.5,2} for the 
left, middle and right plots. The dark shaded lines correspond to non-Markovian walks while 
the white lines indicate non-Markovian jumps to previously visited positions. The bottom 
plot superimposes all of the top three simulations on one graph. 
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Figure 5: The first three subplots in the top row correspond to non-Markovian DODE 
simulations with a = {0.8, 1.3, 1.8} and different values of /3: = .999, = .99 and j3 = .999 
for the left, middle and right plots. The dark shaded lines correspond to non-Markovian 
walks while the white lines indicate non-Markovian jumps to previously visited positions. 
The bottom plot superimposes all of the top three simulations on one graph. 
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